home *** CD-ROM | disk | FTP | other *** search
/ BBS in a Box 3 / BBS in a box - Trilogy III.iso / Files / Prog / U-Z / VideoToolBox Folder / Utilities / Quick3 / TestPsychometricFit.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-02-23  |  4.2 KB  |  130 lines  |  [TEXT/KAHL]
  1. /*
  2. TestPsychometricFit.c
  3. Copyright (c) 1990-1992 Denis G. Pelli
  4. This is a simple driver to show that PsychometricFit.c works.
  5.  
  6. I use the Weibull psychometric function to provide the probability at each contrast
  7. and the Numerical Recipes bnldev() binomial deviates function to simulate the appropriate
  8. number of trials at that contrast.
  9.  
  10. Then I call PsychometricFit() and ask it to fit the Weibull() function to the data, and
  11. then I print out the results. As one would expect, the fits are very good. It seems
  12. reasonable to reject fits at the 5% significance level. So expect to reject about
  13. 5% of your fits even if all the assumptions of this model are correct.
  14.  
  15. Optionally, also compares QUICK3 with QUEST. (Note: QUEST is not part of the 
  16. VideoToolbox. Sorry.)
  17.  
  18. HISTORY:
  19. 4/6/90    dgp    wrote it. Seems to work fine for all cases, 1 to 100,000,000 trials per
  20.             contrast, and few and many contrasts.
  21. 11/18/92 dgp added comparison with QUEST.
  22. 2/20/93    dgp    added call to Require().
  23. */
  24. #include "VideoToolbox.h"
  25. #define QUEST 0
  26. #if QUEST
  27.     #include "Quest.h"
  28. #endif
  29. #include "Quick3.h"
  30. #include <assert.h>
  31. #include <time.h>
  32. #include <nr.h>        /* Numerical Recipes in C */
  33.  
  34. void main(void)
  35. {
  36.     dataRecord data,monotonicData;
  37.     contrastRecord *cPtr;
  38.     paramRecord params;
  39.     int i;
  40.     double chiSquare,weibullLL,monotonicLL;    /* log likelihood */
  41.     int chiSquareDF,weibullDF,monotonicDF;    /* degrees of freedom */
  42.     double p,logC,range;
  43.     int idum;
  44.     int j,cond=0;
  45.     #if QUEST
  46.         Quest *q;
  47.     #endif
  48.     double mode;
  49.     
  50.     Require(0);
  51.     params.logAlpha=0;
  52.     params.beta=3;
  53.     params.gamma=0.111173; 
  54.     params.delta=0.01;
  55.     
  56.     idum=time(NULL);
  57.     printf("Random seed %d\n",idum);    /* So we can reproduce interesting cases */
  58.     data.contrasts=5;
  59.     range=0.5;                            /* log contrast range, centered on logAlpha */
  60.     range=2;
  61.     for(i=0;i<data.contrasts;i++){
  62.         logC=params.logAlpha + range*(i/(data.contrasts-1.0) - 0.5);
  63.         data.c[i].contrast=pow(10.0,logC);
  64.         data.c[i].trials=10;                /* trials at this contrast */
  65.         data.c[i].correct=bnldev(Weibull(data.c[i].contrast,¶ms),data.c[i].trials,&idum);
  66.     }
  67.     #if QUEST
  68.         /* Quest Parameters */
  69.         q=(Quest *)malloc(sizeof(Quest));
  70.         assert(q!=NULL);
  71.         q->nConds=1;    
  72.         q->nLevels=600;
  73.         q->nTrials=0;
  74.         q->grain=0.01;    /* step size of grid, in log contrast */
  75.         q->initialSD=1;
  76.          q->nResponses=2;
  77.          q->quantileOrder=NAN;
  78.          q->fakeIt=0;
  79.          q->function=WeibullPResponse;
  80.         q->beta=params.beta;
  81.         q->gamma=params.gamma;
  82.         q->delta=params.delta;
  83.         q->epsilon=0.0;
  84.         cond=0;
  85.         q->guess[cond]=0;
  86.         QuestOpen(q);
  87.         for(i=0;i<data.contrasts;i++){
  88.             logC=log10(data.c[i].contrast);
  89.             for(j=0;j<data.c[i].correct;j++)QuestUpdate(q,cond,logC,1);
  90.             for(;j<data.c[i].trials;j++)QuestUpdate(q,cond,logC,0);
  91.         }
  92.         q_removePrior(q->qConds[cond]);
  93.         mode = q_mode(q->qConds[cond]);
  94.         QuestClose(q);
  95.     #endif
  96.     printf("Testing the function PsychometricFit.\n");
  97.     printf("Simulating an observer with a Weibull psychometric function.\n");
  98.     weibullDF=1;    /* number of parameters to be adjusted in fitting */
  99.     printf("The simulated data will be fit using %d degrees of freedom.\n",weibullDF);
  100.     printf("Observer: ");
  101.     printf("logAlpha%6.2f, beta%4.1f, gamma%5.2f, delta%5.2f\n",params.logAlpha,params.beta,params.gamma,params.delta);
  102.     p=PsychometricFit(¶ms,&Weibull,&data,&weibullLL,weibullDF,&chiSquare,&chiSquareDF);
  103.     printf("Fit:      ");
  104.     printf("logAlpha%6.2f, beta%4.1f, gamma%5.2f, delta%5.2f, significance%5.2f\n",params.logAlpha,params.beta,params.gamma,params.delta,p);
  105.     #if QUEST
  106.         printf("QUEST mode %.2f\n",mode);
  107.     #endif
  108.  
  109.     /*
  110.     We're done, but just to show off, let's print out everything that anyone
  111.     could possibly want. In real life I would skip this junk.
  112.     */
  113.     monotonicData=data;
  114.     MonotonicFit(&monotonicData,&monotonicLL,&monotonicDF);    /* overwrites data with fit */
  115.     printf("\ncontrast Trials Right  Ratio  Weibull Monotone\n");
  116.     for(i=0;i<data.contrasts;i++){
  117.         cPtr=&data.c[i];
  118.         printf("%6.3f   %5ld %5ld %7.3f %7.3f %7.3f\n",
  119.             cPtr->contrast,cPtr->trials,cPtr->correct,
  120.             cPtr->correct/(double)cPtr->trials,
  121.             Weibull(cPtr->contrast,¶ms),
  122.             monotonicData.c[i].correct/(double)monotonicData.c[i].trials
  123.             );
  124.     }
  125.  
  126.     chiSquare=2.0*(monotonicLL-weibullLL);
  127.     chiSquareDF=monotonicDF-weibullDF;
  128.     p=PChiSquare(chiSquare,chiSquareDF);
  129.     printf("\nChi square %.1f with %d degrees of freedom, yielding a significance of %.2f\n",chiSquare,chiSquareDF,p);
  130. }